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Abstract 



Three-nucleon forces (3NF) are investigated from two-flavor lattice QCD simula- 
tions. We utilize the Nambu-Bethe-Salpeter (NBS) wave function to determine two- 
nucleon forces (2NF) and 3NF in the same framework. As a first exploratory study, we 
extract 3NF in which three nucleons are aligned linearly with an equal spacing. This 
is the simplest geometrical configuration which reduces the huge computational cost 
of calculating the NBS wave function. Quantum numbers of the three-nucleon system 
are chosen to be (J, J p ) = (1/2, l/2 + ) (the triton channel). Lattice QCD simulations 
are performed using Nf = 2 dynamical clover fermion configurations at the lattice 
spacing of a = 0.156 fm on a 16 3 x 32 lattice with a large quark mass corresponding to 
= 1.13 GeV. We find repulsive 3NF at short distance in the triton channel. Several 
sources of systematic errors are also discussed. 
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§1. Introduction 

The nuclear force is an essential input in studying the nuclear many-body problems. In 
particular, the "realistic nuclear potentials" between two nucleons (2N), which can success- 
fully reproduce the vast amount of the 2N scattering data, have been used in ab-initio nuclear 
calculations. It turns out, however, that these two-nucleon forces (2NF) generally underes- 
timate the experimental binding energies of light nucleiP'^ 1 this indicates the necessity of 
three- nucleon forces (3NF). High precision deuteron-proton elastic scattering experiments at 
intermediate energies have also shown a clear indication of 3NFP™ 

There also exist various phenomena in nuclear physics and astrophysics where 3NF may 
play an important role. Some examples include (i) the cross section of the backward scatter- 
ing angles in nucleus-nucleus elastic scattering,-' (ii) the anomaly in the oxygen isotopes near 
the neutron drip-line,® and (iii) the nuclear equation of state (EoS) at high density relevant 
to the physics of neutron stars P Universal short-range repulsion for three baryons (nucleons 
and hyperons) has also been suggested in relation to the maximum mass of neutron stars 
with hyperon coreP® 

Despite of its phenomenological importance, microscopic understanding of 3NF is still 
limited. Pioneered by Fujita and Miyazawap^ 1 the long range part of 3NF has been mod- 
eled by the two-pion exchange (27rE)p^ particularly with the Z\-resonance excitation. This 
27rE-3NF component is known to have an attractive nature at long distance. An additional 
repulsive component of 3NF at short distance is often introduced in a purely phenomenolog- 
ical wayJ-^ 1 An approach based on the chiral effective field theory (EFT) is quite useful to 
classify and parametrize the two-, three- and more-nucleon forces .QSJiUDilS While unknown 
low-energy constants cannot be determined within its framework but are obtained only by 
the fitting to the experimental data, the obtained parameters are used in nuclear physics 
calculations.^ 1 

To go beyond phenomenology, it is most desirable to determine 3NF directly from the 
fundamental degrees of freedom (DoF), the quarks and the gluons, on the basis of quantum 
chromodynamics (QCD). In this letter, we make a first exploratory study of first-principle 
lattice QCD calculation of 3NFn 

As for the 2NF from lattice QCD, an approach based on the Nambu-Bethe-Salpeter 
(NBS) wave function has been proposedP^'^ 1 Resultant (parity-even) 2NF in this approach 
are found to have attractive wells at long and medium distances and central repulsive cores 
at short distance. The method has been extended also to the hyperon-nucleon (YN) and 



*' We note here that a completely different approach based on holographic QCD leads to a nuclear 
matrix model whose solution gives repulsive 3NF at short distance.^ 
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hyperon-hyperon (YY) interactions.^'^'^'GS From the numerical point of view, the ex- 
tension of the above methodology to three-nucleon (3N) systems is rather challenging due 
to huge number of DoF of the NBS wave function for 3N systems. Also, we need to identify 
genuine 3NF by subtracting out the 2NF in both parity-even channel and parity-odd channel. 
The latter channel has not yet been extracted in the present lattice QCD simulations. 

In this paper, as a first exploratory study, we consider the triton channel, / = 1/2, 
J p = l/2 + , as the quantum numbers of the 3N system. In this case, we can determine 
3NF even without the knowledge of parity-odd 2NF. In order to reduce the computational 
cost, we study 3NF through a configuration of three nucleons aligned linearly with an equal 
spacing, which is the simplest geometrical configuration (a linear setup) for the 3N system. 

Lattice QCD simulations are carried out with Nf = 2 dynamical clover fermion con- 
figurations. Preliminary accounts of this study have been given in Refs. I24"|) . l25]) . It is in 
order here to mention that lattice QCD simulations for three- and four- baryon systems have 
been also performed in Refs. l2"6j) . [2"Tj) . They, however, focus on the energies of the multi- 
baryon systems, and extracting 3NF, which could be used in future ab-initio nuclear physics 
calculations, is currently beyond their scope. 

This paper is organized as follows. In Section [21 we describe the basic formulation to 
study 3N systems, as well as its challenges. In Section [3j we develop a framework to obtain 
3NF in the triton channel using only parity-even 2NF. The advantage of the "linear setup" 
for the three-dimensional (3D) coordinate configuration is also discussed. In Section HJ we 
explain our lattice QCD setup, and numerical results of 3NF are presented in Section [51 
Section [6] is devoted to summary and concluding remarks. In Appendix [A], we present a 
complementary study using effective 2NF in 3N systems. 

§2. Formulation for three-nucleon systems 

To begin with, we briefly review the framework for the calculation of 2N potentials.^ 1 
We consider the following effective Schrodinger equation, 

- ^m(f) + J dfU 2N (ry)iP 2N (f) = E 2N ^j 2N (r), (2-1) 

where \i = m N /2 and E 2N are the reduced mass and ground-state energy of the 2N system 
in the center-of-mass frame, respectively, and U 2 n(t, f) is the non-local energy- independent 
2N potential. The equal-time NBS wave function ip 2 N{r) is extracted from the four-point 
correlator as 

G 2N (r,t-t ) = ^^(0\(N(R + ^N{R)){t)W^(to)\0), (2-2) 

R 
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-* A 2N ij 2N e- E ^ t - t0 \ A 2N = (E 2N \(N'N')\0), (2-3) 



i, 2N {r) = (0\N(f)N(0)\E 2N ), (24) 

where \E 2 ^) denotes the ground state vector of the 2N system, N (JV') the nucleon operator 
in the sink (source), and spinor/flavor indices are suppressed for simplicity. In the practical 
calculation, we perform the derivative expansion for the non-locality of the potential,^ 

UMr,?) = [V c (r) + V T (r)S 12 + V LS (r)L ■ S + (9(V 2 )] 8(f - 0, (2-5) 

where Vc, Vp and Vls are the central, tensor and spin-orbit potentials, respectively In 
Ref. 129]) . the validity of this expansion is examined, and it is shown that the leading terms, 
Vc and Vr, dominate the potential at low energies. 

The extension to 3N systems is described as follows. We consider the NBS wave function 
^zn (r, P) extracted from the six-point correlator as 

G 3N (r,p,t-t Q ) = ^2 l {Q\{N{xx)N{x 2 )N{x 3 )) (t) JWN^Wj(t )\0), (2-6) 

R 

-> A 3N ij 3N (r,P)e- E * N{t - t0 \ A 3N = (E m \(N> N> N>)\0) , (2-7) 



t>t 

^ 3 7v(r,p) = (OIN^N^N^IEsn), (2-8) 

where E 3N and \E 3N ) denote the energy and the state vector of the 3N ground state, re- 
spectively, and R = [x\ + x 2 + x 3 )/3, r = x\ — x 2 , p = x 3 — (x x + x 2 )/2 the Jacobi 
coordinates. Assuming that the NBS wave function has a proper asymptotic behavior at 
outer non-interacting regions, as was proven in the case of 2N systems,^'^'^'^'^ 1 we 
consider the following Schrodinger equation of the 3N system with the derivative expansion 
of the potentials, 



l -W 2 r - -^Vj + v Mn 3 ) + V 3NF (f, p) 



i>w{r, p) = E 3N ip m (r, p), (2-9) 



where V 2 jsr(fij) with fij = x j — Xj denotes the 2NF between (z,j)-pair, V 3 Np(f,p) the 3NF, 
/i r = m N /2, Hp = 2m N /3 the reduced masses. 

3NF can be determined as follows. We first calculate iP 3 n(^,P) for all r, p and obtain 
the total potential of the 3N system through Eq. fl2-9p . We also perform (separate) lattice 
simulations for genuine 2N systems and obtain all necessary V 2 N(fij). We then extract 
V 3 nf{t.,p) by subtracting ^^^(r^) from the total potential. The extension to four- and 
more-nucleon forces can be immediately understood. Note that potentials determined in this 
way reproduce the energy of the system by construction. 
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Important remark is that 3NF are always determined in combination with 2NF, and 3NF 
alone do not make too much sense. The comparison between lattice 3NF and phenomeno- 
logical 3NF can be done only at a qualitative level. Rather, our purpose is to determine 
two-, three-, (more-) nucleon forces systematically, and provide them as a consistent set. 

In the practical calculation with the above procedure, the computational cost is enormous: 
Note that adding one more nucleon to the system corresponds to adding three (valence) 
quarks, each of which has three color and four Dirac spinor DoF. In addition, the number of 
diagrams to be calculated in the Wick contraction tends to diverge with a factor of N u \ x iVJ, 
where N u (Nd) are numbers of up (down) quarks in the system. We develop several techniques 
to reduce these computational costs. For instance, we take advantage of symmetries (such 
as isospin symmetry) to reduce the number of the Wick contractions. We also employ 
the non-relativistic limit for the nucleon operator in the source to reduce the number of 
spinor DoF. Similar techniques are (independently) developed in Ref . . Nevertheless, the 
simulation cost remains quite enormous, compared to the current computational resources, 
so we explore an efficient way to restrict the spacial DoF of r and p, as will be explained in 
the next section. 

Even after we manage to calculate ^ 3 jv(r*, p), the identification of genuine 3NF remains 
nontrivial. In general, in order to subtract the contributions from V^at in Eq. (12-91) . we 
need V2N m both of parity-even and parity-odd channels, since a 2N-pair inside the 3N 
system could be either of positive or negative parity. In practice, however, the parity- 
odd V2N have not been obtained yet in lattice QCD. Note that the partial wave expansion 
cannot be performed here, since we calculate only restricted DoF of r and p to reduce the 
computational cost. In the next section, we establish a framework to extract 3NF without 
referring to parity-odd 2NF. 

§3. Identification of three-nucleon forces and the linear setup for 

three-nucleon locations 

We study the 3N system in the triton channel, and consider Eq. (I2-9|) with fixed 3D 
coordinate configurations of 3N. In particular, we choose small |rj, |p| configurations. This 
corresponds to locating 3N close to each other, and the 3NF effect is expected to be enhanced. 

To proceed, the lack of the information of parity-odd 2NF on the lattice prevents us from 
identifying 3NF, as was explained in Section [2j In order to resolve this issue, we consider 
the following channel, 



(3-1) 
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which is anti-symmetric in spin/isospin spaces for any 2N-pair. Combined with the Pauli- 
principle, it is automatically guaranteed that any 2N-pair couples with even parity only. 
Therefore, parity-odd 2NF vanish in (tps\H\ip3N)i where H is the Hamiltonian of the 3N 
system, and we can extract 3NF unambiguously using only parity-even 2NF. Note that no 
assumption on the choice of 3D-configuration of r, p is imposed in this argument, and we 
can take advantage of this feature for 3NF calculations with various 3D-configuration setup, 
which will be our subject in future simulations. 

As noted in Section [21 it is important to restrict the spacial DoF of r and p efficiently, 
in order to reduce the computational cost of ^^{r,p). In this paper, we propose to use the 
linear setup with p = 0. In this setup, three nucleons are aligned linearly with equal spacings 
of T2 = \r2\ with r*2 = r/2. 

The advantage of the linear setup is understood as follows. Because of p = 0, the third 
nucleon is attached to (1, 2)-nucleon pair with only S-wave. Considering that the total 
3N quantum numbers are / = 1/2, J p = l/2 + (triton channel), the wave function can be 
completely spanned by only three bases, labeled by the quantum numbers of (l,2)-pair as 
2S+1 Lj = 1 Sq, 3 Si, 3 -Di- Therefore, the Schrodinger equation can be simplified to the 3x3 
coupled channel equations with the bases of ipis , i> 3 Si, V^zv The reduction of the dimension 
of bases is expected to improve the S/N as well. It is worth mentioning that considering 
the linear setup is not an approximation: Among various geometric components of the wave 
function of the ground state, we calculate the (exact) linear setup component as a convenient 
choice to study 3NF. While we can access only a part of 3NF from it, we plan to extend the 
calculation to more general geometries step by step, toward the complete determination of 
the full 3NF. 

We examine the explicit form of the potential matrix of 2NF. At the leading order of the 
derivative expansion, 2NF is written in terms of central potentials Vq S and tensor potentials 
Vf s with isospin I and spin S, {Vg°, Vg°, Vg 1 , V^, V$\ V^ 1 }, where the label "2iV" is 
omitted for simplicity. An explicit calculation gives us 
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i<j 
( 



p=0 



V™(r 2 ) + VS\r 2 ) 
+ \V™{r) + \V«\r) 


|^°(r 2 ) - \V£\r 2 ) 
-\Vgir) + \VS\r) 


> 

-2VjP(r 2 ) + 2V$ 1 (r) 




\V™{r 2 ) + \V™{r 2 ) + \VS\r 2 ) 
+lV?{r 2 ) + \V™{r) + \V$\r) 


VSHr 2 ) - 3F r u (r 2 ) + 2V T 01 (r) 




V°\r 2 )-W^{r 2 ) + 2V° 1 {r) 


-3y T u (r 2 ) + Vg x (r) - 2V^(r) 



\ 



(3-2) 



where we span the spaces with the rotated bases given by (ips, i>M , i j3 D 1 ) T , where ips hi 
Eq. ([34]) is shown to be tp s = ^(-V^So + and tp M = ■^(+i' 1 s Q + tpwj. Note 

that neither of parity-odd 2N potentials, V® , Vq , V^ 1 , appear in the first row/column in 
Eq. (13 -2p . which is a consequence of the discussion above. 

§4. Setup of lattice QCD simulations 

We employ Nf = 2 dynamical configurations with mean field improved clover fermion 
and renormalization-group improved gauge action generated by CP-PACS Collaboration.® 
We use 598 configurations at (3 = 6/g 2 = 1.95, where the lattice spacing of a -1 = 1.269(14) 
GeV is determined from the rho meson mass,® and the lattice size of V = L 3 x T = 16 3 x 32 
corresponds to (2.5 fm) 3 box in physical spacial size. For u, d quark masses, we take the 
hopping parameter at the unitary point as k u = = k u< i = 0.13750. The masses of pion, 
nucleon and Delta are obtained as m^a = 0.8934(4), tunci = 1.694(1) and m^a = 1.820(2) 
respectively, corresponding to m n = 1.13 GeV, m N = 2.15 GeV, = 2.31 GeV. 

We use the wall quark source with Coulomb gauge fixing, and periodic (Dirichlet) bound- 
ary condition is imposed in spacial (temporal) direction. Due to the enormous computational 
cost, we can perform the simulations only at a few sink time slices. Looking for the range 
of sink time where the ground state saturation is achieved, we carry out preparatory simu- 
lations for effective 2NF in the 3N system in the triton channel. We find that effective 2NF 
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at different sink time slices agree with each other as long as (t — to) /a > 7. Being on the 
safer side, we decide to perform linear setup calculations at (t — t )/a = 8 and 9. For details 
of the study of effective 2NF, see Appendix [S] 

In order to enhance the statistics, we perform the measurement at 32 source time slices 
for each configuration, and the forward and backward propagations are averaged. The results 
from both of total angular momentum J z = ±1/2 are averaged as well. The error is estimated 
by the jackknife method with the bin size of 13, and the auto-correlation is found to be small 
by comparing the results with those from different bin sizes. We carry out the simulation at 
eleven physical points of the distance r 2 = r/2 with the linear setup. More specifically, we 
take f 2 /a= (0,0,0), (1,0,0), (1,1,0), (1,1,1), (2,0,0), (1,1,2), (2,2,0), (3,0,0), (2,2,2), 
(3,3,0), (3,3,3), and additional lattice points obtained by cubic rotation from the above. 
These points are chosen so that they can make a good sampling of r 2 as a whole, and the 
number of lattice points in cubic rotation is small. Furthermore, in order to suppress the 
finite volume artifact, points with relatively large r 2 (r 2 ^ 0.5 fm) are located in off-axis 
direction. 

In the calculation of the four- and six-point correlators (Eqs. (12-21) and (I2-6P ). we have a 
freedom to choose any nucleon interpolating operator in the sink and source, in so far as it 
(strongly) couples to a nucleon state. Since a NBS wave function is defined through a sink 
operator N, a different potential is obtained from a different operator N. Note, however, 
that physical observables calculated from these different potentials, such as phase shifts 
and binding energies, are unique by construction.^ In this sense, one can consider that 
choosing N corresponds to choosing the "scheme" to define the potentialP^'^ We choose 
N so that the non- locality of the obtained potential is small, achieving a reasonable S/N for 
the potential at the same time. To this end, we have found that the standard nucleon (sink) 
operator iV = N std = e a b c (qaC , y 5 qb)qc satisfies these criteria in the 2N studyP^ 1 We therefore 
employ the same sink operator in the 3N study as well, so that 2NF and 3NF are determined 
on the same footing. We still have an additional freedom to choose a source nucleon operator 
N'. Since this choice does not affect the NBS wave functions nor the potentials, our main 
concern here is the reduction of the computational cost. In the 3N study, we employ the 
non-relativistic limit operator, N' = N nr = € a bc(qaCl5Pnrqb)PnrQc with P nr = (1 + 74)/2, 
so that the spinor DoF is reduced and the faster computation is achieved. In the 2N case, 
we perform explicit lattice simulations for both of N' = N nr choice and N' = iV std choice. 
Theoretically, these potentials should agree at t — > oo, and we observe that these are actually 
consistent within statistical fluctuations at our time region of (t — to)/a = 8 and 9. Hereafter, 
we employ the source operator of N nr in the 2N study as well. When extracting 3NF, we 
expect that this choice leads to a better cancellation between 2NF and 3NF for the statistical 
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Fig. 1. Parity-even 2NF obtained at (t — to)/a = 8. V^ denotes the central force in 1 Sq channel, 
and Vq and V^ 1 the central and tensor forces in 3 S\- 3 Di channel, respectively. 



fluctuations as well as systematic error by the excited state contaminations. 

§5. Numerical results of lattice QCD simulations 

We first show numerical results for 2N systems. Shown in Fig. [T] are the parity-even 2NF 
of Vq°, Vq 1 and V? 1 obtained at (t — to) /a = 8, where we have checked that the sink time 
dependence is small as long as (t — to)/a > 7. Vq is obtained from the single channel analysis 
of the state, while Vg 1 and Yg 1 are obtained from the 2x2 coupled channel analysis 
in the 3 S\- 3 Di state. Each potential consists of the kinetic term (Laplacian term) and the 
energy term in Eq. (12- 1|) . The latter is estimated to be ~ —3 MeV and E®^ ~ —3 MeV 
in 1 So and 3 S±- 3 D 1 channels, respectively, using the condition that the potential approaches 
zero at large r (the method called "from V" in Ref. I3"2l ). 

We observe central repulsive cores at short distance, r < 0.2-0.3 fm, and attractive wells 
at middle through long distances, 0.3 ^ r <, 0.8 fm. The qualitative features are same as 
our previous lattice simulations,^'^ 1 while quantitative strength and range of potentials 
are associated with the particular lattice setup in this calculation, e.g., quark masses and a 
lattice spacing. 

Let us turn to the results for the 3N system in the triton channel. Each of ipsi ^m, ^ 3 Di 
is extracted by the projection in spin and flavor spaces. Simultaneously, the projection to 
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Fig. 2. 3N wave functions at (t — to) /a = 8. Circle (red), triangle (blue), square (green) points 
denote ips, ipMt ^s^, respectively. r 2 = r/2 (with offset for visibility) is the distance between 
the center and edge in the linear setup. 



the representation of the cubic group (which corresponds to the S-wave projection in the 
continuum limit) is performed on if)$ an d i^m, while the subtraction of the Af representation 
is performed on ip3£ >1 . In ip3£ >1 , there are various Yj=2,m spherical harmonics components. We 
find that they are generally consistent with each other by making a division with a proper 
Y 2jm . This indicates that the contamination from higher waves in ip3 Dl are negligible, and 
that the artifact due to the breaking of the rotational symmetry to the cubic symmetry is 
small O 

In Fig. [21 we plot the radial part of each wave function of ij)g, ipM, V^Di obtained at 
(t — to) /a = 8. Here, we normalize the wave functions by the center value of ips{ r 2 — 0). 
What is noteworthy in Fig. [2] is that the wave functions are obtained with good precision. 
This is quite nontrivial, since the S/N generally gets worse when the system of concern is 
composed of more quarksP^ 1 Furthermore, by fixing the 3D-configuration of 3N, each data 
point is obtained from effectively fewer number of statistical samplings. This is in contrast 
to, e.g., calculations of the energy, where there exists an effective gain on number of statistics 
by summations over the lattice volume for each nucleon. 

*) Small discrepancy is observed only at r2/a — y2j which signals a splitting between E and 
representations. This artifact, however, makes negligible effect on results of 3NF. 
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Fig. 3. The effective scalar-isoscalar 3NF in the triton channel with the linear setup obtained at 
(t-t )/a = 8. 

In Fig. [2J we observe that tps overwhelms the wave function, and ipM, ip^Dx are much 
smaller by one to two orders of magnitude. This indicates that higher partial wave com- 
ponents in ips are also strongly suppressed, and the wave function is completely dominated 
by the component with which all three nucleons are in S-wave. We consider that this is 
mainly originated by the large quark mass employed in this calculation. For instance, 1^3^ 
component is expected to be suppressed due to smaller tensor forces with heavier pions. 
In fact, in the case of spin-triplet 2N system, lattice simulations with several quark masses 
explicitly show that the D-wave component is smaller at larger quark masses.^ 1 

We determine 3NF by subtracting V2N from the total potentials in the 3N system. As 
is explicitly shown in Eq. (13-21) . we have only one channel, (ips\H\ip3N) , which is free from 
the parity-odd V27V. Correspondingly, we can determine only the sum of the 3NF matrix 
elements, (^s\V3nf\^3n) = (^sl^Wl^s) + {^s^nf^m) + (i>s\VwF\i>a Dl ), or one type 
of spin/isospin functional form for 3NF. In this paper, 3NF are effectively represented in a 
scalar-isoscalar functional form. This form is an efficient representation, since ips overwhelms 
the wave function and thus K^sI^WfItAs)! ~> \(^Ps\V3Nf\^Pm)\, \{^s\V3Nf\^ Di )\ is expected. 
Note also a scalar-isoscalar functional form is often employed for the short-range part of 3NF 
in phenomenological models.^ 1 

When examining 3NF, we have to include r 2 -independent shift by energies, in addition 
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Fig. 4. 3NF as explained in Fig. [3J Circle (red) points are obtained at (t — to)/a = 8, and triangle 
(blue) points are obtained at (t — to) /a = 9. 



to kinematic terms (Laplacian part in Eqs. (j2-ip and (I2-9P ). As was done for the 2NF, we 
employ the "from V" method,^ and obtain E%n — —3 MeV from the long-range behavior 
of the effective 2NF in the 3N system (see Appendix |A] for details.) Combined with the 
energies of 2N, the total energy shift for the 3NF amounts to 5e = E 3N — ?>/2E\° N — ?>/2E^ N 
~ 5 MeV. Note that the precise determination of energies is still a challenging issue, and 
we estimate that 5e suffers from < 10 MeV systematic error. This uncertainty, however, 
does not affect the following discussions much, since Se merely serves as an overall offset. In 
Fig. [31 we plot the results for the effective scalar-isoscalar 3NF at (t — t )/a = 8. In order to 
check the dependence on the sink time slice, we compare 3NF from (t — to) /a = 8 and 9 in 
Fig. HI While the results with (t — to) /a = 9 suffer from quite large errors, they agree with 
each other within statistical fluctuations. 

Fig. [3] shows that 3NF are small at the long distance region of T2- This is in accordance 
with the suppression of 27rE-3NF by the heavy pion. At the short distance region, on the 
other hand, we observe an indication of repulsive 3NF. Note that a repulsive short-range 
3NF component is phenomenologically required to explain the properties of high density 
matter. Since multi- meson exchanges are strongly suppressed by the large quark mass on 
a lattice, the origin of this short-range 3NF may be attributed to the quark and gluon 
dynamics directly. In fact, we recall that the short-range repulsive (or attractive) cores in 
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Fig. 5. 3NF as explained in Fig. [3l Circle (red) points are obtained from Vg td , and triangle (blue) 
points from V? . See the text for details. 



the generalized baryon-baryon potentials are systematically calculated in lattice QCD in the 
flavor SU(3) limit and the lattice results are found to be well explained from the viewpoint 
of the Pauli exclusion principle in the quark levelJ^'^ 1 '® i n this context, it is intuitive to 
expect that the 3N system is subject to extra Pauli repulsion effect, which could be an origin 
of the observed short-range repulsive 3NF. Further investigation along this line is certainly 
an interesting subject in future. 

Discussions on several systematic errors are in order. First, one may worry about the 
discretization error, since the nontrivial results are obtained in the short-range part of 3NF. 
In particular, the kinetic terms (Laplacian part in Eqs. ( 12 -ip and (12-91) ) could suffer from a 
substantial effect, since they are calculated by the finite difference as 

V 2 /(^) = V s 2 td /(:r) = l£ ^ x + a *) + K x ~ a *) " 2 ^ ■ t 5 ' 1 ) 

cl 

i 

In order to estimate this artifact, we carry out additional analyses employing the improved 
Laplacian operator for both of 2N and 3N systems, 

i 

+16(f(x + a t ) + f(x - a,)) - 30/(x))] . (5-2) 
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Shown in Fig. |5] is the comparison between the results of 3NF from Vg td and Vf mp at (t — 
to) /a = 8. We observe that they are consistent with each other, and conclude that the 
discretization artifact of 3NF in Laplacian operator is small. We, however, remark that this 
study probes only a part of discretization errors, and an explicit calculation with a finer 
lattice is desirable. Actually, the analysis with operator product expansion^ shows that 
2NF of V2jv(r) tend to diverge as r — > 0, so significant discretization artifact is expected 
around r = 0. A simulation with a finer lattice is currently underway, and we plan to report 
it in future publications. 

Second, the finite volume artifact could be substantial, considering that three nucleons 
are accommodated in (2.5 fm) 3 spacial lattice box. Such artifact, however, would mainly 
appear in large r2 region, and we avoid it as much as possible by focusing on the short- 
range part of 3NF. Furthermore, recall that points with relatively large r2 (V2 0.5 fm) are 
carefully chosen so that they are located in off-axis directions. We also note that the size of 
(2.5 fm) 3 is quite large for the heavy pion in this calculation, namely, m n L is as large as 14. 
Of course, a quantitative estimate requires calculations with larger volumes, which we defer 
to future studies. 

Third, we study the 3NF at leading order in the derivative expansion. Although it has 
been shown that higher order terms are small in (parity-even) 2NFJ20 1 tiny modifications in 
2NF could make a nonnegligible effect on 3NF. On this point, we recall that Fig. [2] indicates 
that the wave function is completely dominated by the component with which all three 
nucleons are in S-wave. Therefore, contaminations from spin-orbit forces, which is the next 
leading order in the derivative expansion, are expected to be small in this lattice setup. It 
would be interesting to see how much higher derivative terms affect 3NF in physically lighter 
quark mass, where, e.g., D-wave components get larger by larger tensor forces. 

Since the lattice simulations are carried out only at single large quark mass point, quark 
mass dependence of 3NF is certainly an important issue to be investigated. In the case of 
2N or generalized baryon-baryon potentials, short-range cores have the enhanced strength 
and broaden range by decreasing the quark mass JI3>[23>[23>[I3 We, therefore, would expect a 
significant quark mass dependence exist in short-range 3NF as well. In addition, long-range 
27rE-3NF will emerge at physically light quark masses. As is well known, such calcula- 
tion poses a significant challenge to lattice QCD, since S/N is quickly ruined by pursuing 
lighter quark mass direction. When one considers systems consisting of more nucleons, more 
severely this issue emerges.® Theoretical developments are highly required, and some of 
our achievements are being reportedP^^ 
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§6. Summary 



We have studied three-nucleon forces (3NF) in the triton channel, I = 1/2, J p = l/2 + , 
in lattice QCD, using Nf = 2 dynamical clover fermion configurations with the size of 
16 3 x 32, the lattice spacing of a = 0.156 fm and a large quark mass corresponding to 
m n = 1.13 GeV. We have utilized the Nambu-Bethe-Salpeter wave function to determine 
two-nucleon forces (2NF) and 3NF on the same footing. We have developed a framework to 
identify 3NF unambiguously, without referring to the information of parity-odd 2NF. When 
extracting 3NF, we have fixed the three-dimensional (3D) coordinate configuration of three 
nucleons. We have proposed to employ the "linear setup" for the 3D-configuration, where 
three nucleons are aligned linearly with equal spacings. With this setup, we have shown 
that the Schrodinger equation can be simplified to the 3x3 coupled channel equations. 
Performing lattice simulations with the linear setup, we have found repulsive 3NF at short 
distance. Several sources of systematic errors have been discussed, and prospects for future 
lattice simulations have been given as well. In particular, it is important to pursue the 
study with lighter quark masses, considering broad impacts of 3NF on various phenomena in 
nuclear physics and astrophysics, e.g., the natures of high density matter and neutron stars. 
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Appendix A 

Effective two-nucleon forces in three-nucleon systems 

We investigate the Schrodinger equation for the 3N system, Eq. f!2-9p . from a viewpoint 
of effective 2NF. More specifically, we take the summation over the location of a spectator 
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Fig. 6. Same as Fig. except that the effec- 
tive 2NF in the triton channel are plotted, 
instead of the genuine 2NF. 
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nucleon N{xz), 



(a-i; 



and define the effective potential between N(xi) and N(x 2 ) through the following effective 
Schrddinger equation, 



(A-2) 



Since the DoF of p is integrated out beforehand, computational cost is reduced drastically. 
This summation is also expected to enhance the S/N. Due to these features, effective 2NF 
serve as a suitable framework for preceding simulations to much more expensive linear setup 
simulations (or any other fixed 3D-configuration simulations.) It can also determine E$n, 
which provides an ^-independent shift of 3NF in the linear setup simulations. 

As the 3N system, we consider the triton channel, I = 1/2, J p = l/2 + . Because the 
spectator nucleon is projected to S-wave, possible quantum numbers between the (effective) 
2N are only 2S+1 Lj = 1 S , 3 S U 3 D 1 . We calculate the effective 2NF Vjg, i.e., the central 
Vq° sS , Vq\ s and the tensor V^\ s potentials, where I and S denote the isospin and spin, 
respectively. 

Lattice simulations are carried out with the same setup as in Section m except for that 
the measurement is taken at 16 source time slices for each configuration, and with wider 
sink time range of 2 < (t — to) /a < 11. In Fig. [BJ we show results for V^{r) obtained at 
(t — to) /a = 8. The results agree with those from different t as long as (t — to)/ a > 7. The 
r-independent shift by the energy in Eq. (IA-2j) is determined from the long range behavior 
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of the potential (the method called "from V" in Ref. 32)). We obtain E^n ~ —3 MeV from 
both of V3e ff and Vq\ s , consistently. Note that the tensor potential does not suffer from 
this overall shift 

It would be of interest to compare effective 2NF V e s(r) and genuine 2NF V2n{ t ) and 
extract the effect of the 3N system. Some of the effects are attributed to genuine 2NF with 
nontrivial 3N correlations, and the others are originated by genuine 3NF. In this way, we 
can indirectly access the effect of 3NF. 

In Fig. [?l we plot V e s(r) — V2n(j) for the tensor potential, V^ 1 . One can see that the 
difference is consistent with zero, and there is no indication of the 3NF effect. The results 
of Kff( r ) — ^2Af( r ) for the central potentials Vq° and V® 1 are similarly consistent with zero, 
with larger statistical errors than the case of V® 1 by a factor of 2-4. 

This (non-)observation can be explained as follows. First of all, the 27rE-3NF component 
is expected to be suppressed with m w = 1.13 GeV in this calculation. Furthermore, A- 
excitation is suppressed by the S-wave projection on a spectator nucleon. We still have 
short-range repulsive 3NF observed in Section It turns out, however, that such 3NF effect 
is obscured due to the volume factor in the summation of Eq. (lA-ip . For instance, if we 
consider effective 2NF with the distance of r ~ 0.2 fm, a spectator nucleon is involved in 
3NF only when it resides within < 0.4 fm from the center of the 2N of concern. Therefore, 
the relative volume factor amounts to <, |7r(0.4fm) 3 /(2.5fm) 3 = 0.017 and thus the 3NF 
affect effective 2NF only <C 1 MeV. This magnitude is comparable to the size of the statistical 
errors in the effective 2N study. 

These discussions in turn show that it is essential to perform the simulation with fixed 
3D-configuration of 3N in order to extract 3NF. 
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